Energy level statistics of interacting trapped bosons 



Barnali Chakrabarti 1,2 , Anindya Biswas 3 . V. K. B. Kota 4 , Kamalika Roy 2 , Sudip Kumar Haldar 2 
Instituto de Fisica, Universidade de Sao Paulo, CP 66318, 05315-970, Sao Paulo, SP Brazil. 
2 Department of Physics, Lady Brabourne College, 
Pl/2 Suhrawardi Avenue, Kolkata 700017, India. 
3 Department of Physics, University of Calcutta, 92 A.P.C. Road, Kolkata 700009, India. 
4 Physical Research Laboratory, Navarangpura, Ahmedabad 380009, India. 

It is an well established fact that statistical properties of energy level spectra are the most efficient 
tool to characterize nonintegrable quantum systems. Statistical behavior of different systems like, 
complex atoms, atomic nuclei, two-dimensional Hamiltonians, quantum billiards and non-interacting 
many bosons have been studied. The study of statistical properties and spectral fluctuation in the 
interacting many boson systems have developed a new interest in this direction. Specially we are 
interested in the weakly interacting trapped bosons in the context of Bose-Einstein condensation 
(BEC) as the energy spectrum shows a transition from the collective to single particle nature with 
the increase in the number of levels. However this has received less attention as it is believed that 
the system may exhibit Poisson like fluctuations due to the existence of external harmonic trap. 
Here we compute numerically the energy levels of the zero-temperature many-boson systems which 
are weakly interacting through the van der Waals potential and are in the 3D confined harmonic 
potential. We study the nearest neighbour spacing distribution and the spectral rigidity by unfolding 
the spectrum. It is found that increase in number of energy levels for repulsive BEC induces a 
transition from a Wigner like form displaying level repulsion to Poisson distribution for P(s). It 
does not follow the GOE prediction. For repulsive interaction, the lower levels are correlated and 
manifest level repulsion. For intermediate levels P(s) shows mixed statistics which clearly signifies 
the existence of two energy scales : external trap and interatomic interaction, whereas for very high 
levels the trapping potential dominates, genarating Poisson distribution. Comparison with mean- 
field results for lower levels are also presented. For attractive BEC near the critical point we observe 
the Shnirelman like peak near s = which signifies the presence of large number of quasi-degenerate 
states. 

PACS numbers: 03.75.Hh, 31.15.Ja, 05.45.Mt, 05.45.Pq 



I. INTRODUCTION 

Although there is no precise definition of the qun- 
tum chaos, but the statistical properties of the energy 
level spectra is often used to characterize level fluctua- 
tion in quantum systems. It is an well established fact 
that for the classically integrable systems the energy lev- 
els are uncorrelated and nearest neighbor level spacing 
distribution follows Poisson statistics [l[ , whereas classi- 
cally chaotic systems are associated with spectral fluctu- 
ations and strong level repulsion between energy levels is 
described in random matrix theory [2-4]. It follows GOE 
(Gaussian orthogonal ensemble) or GUE (Gaussian uni- 
tary ensemble) of random matrices depending whether 
the Hamiltonian has time-reversal symmetry or not @, 0] . 
The spectral properties of many different many fermion 
quantum systems like atoms and atomic nuclei and also 
quantum billiards have already been studied [3MTT| . P{s) 
distribution of nuclear data ensemble agrees very well 
with GOE and in the atomic spectra the nearest-neighbor 
spacing distribution shows Wigner type. In addition, re- 
cently there are some studies of spectral properties of 
non-interacting many particle (fermions or bosons) sys- 
tems [IH and interacting boson systems [l(| EH, E34I3 • 
In the present work we undertake to study the 
quantum mechanical spectra, the statistical behavior of 
weakly interacting many-boson systems with an external 



harmonic confinement. This study is specially interest- 
ing for several reasons. Firstly it directly corresponds 
to the Bose-Einstein condensation in dilute atomic vapor 
0, E3 ■ Secondly due to the presence of external con- 
finement the energy spectrum shows a transition from 
collective to single particle nature [HI, EH . Apparently it 
appears that the system will exhibit most common Pois- 
son statistics of integrable systems as at the near zero 
temperature the interaction energy is small compared to 
the trap energy. However from the earlier studies of dif- 
ferent statistics, thermodynamic and dynamic properties 
of the system it is an established fact that interatomic 
interaction plays an imp ortant role even in the weakly 
interacting Bose gas [16l - [F9l ]. Naturally it leads us to be 
more curious in the study of level spacing distribution. 
It seems to contradict the usual expectation based on 
RMT. This is a different type of system where two energy 
scales coexist. One is the external trap which is charac- 
terized by the trap energy Hlo (to is the external trap 
frequency). The other one is the interatomic interaction 
which is characterized by Na s ; where N is the number of 
bosons in the trap and the properties of zero-temperature 
BEC are essentially characterized by the s-wave scatter- 
ing length a s . Thus the study of spectral statistics of 
such a realistic system may provide exciting information 
on the level correlation and may disagree the universal 
hypothesis of Bohigas, Giannoni and Schmit Due to 
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existence of two energy scales, the system does neither 
obey the regular Poisson distribution nor the GOE for a 
strongly chaotic system. We observe new features from 
the following study. Quantum mechanical spectra un- 
dergoes a transition, as a function of number of energy 
levels. It has already been observed both experimentally 
and theoretically that low-lying levels are strongly influ- 
enced by interatomic interaction [lU, [l9[ . These levels are 
highly correlated and we observe close to Wigner distri- 
bution. The intermediate levels show a mixed statistics 
which is the overlap of both Poisson and Wigner distribu- 
tion which clearly signify the co-existence of two energy 
scales. Thus the choice in the number of levels has a great 
influence in the statistical properties. For higher levels 
(much above the chemical potential) the energy spectra 
is strictly dominated by the harmonic confinement and 
energy levels are almost equidistant by the amount Huj, 
similar to the non interacting harmonic oscillator and it 
generates Poisson distribution. To the best of our knowl- 
edge there is neither any systematic calculations nor any 
rigorous derivation in this direction. Here we tackle the 
problem by solving the trapped many-body system b y an 
ab initio but approximate many-body technique [20N22| . 
As it is complicated many-body problem, it is hard to 
present analytic studies. However our numerical study is 
also important to investigate statistical behavior of such 
a realistic condensate. 

The paper is organized as follows. Sec II deals with 
the many-body technique which basically calculates the 
many-body effective potential. Choice of interaction and 
detailed calculation of energy levels are presented in Sec 
III. Sec IV deals with several statistical tools and numer- 
ical results. Sec.V concludes the summary. 



II. MANY-BODY TECHNIQUE 

We start with the Hamiltonian of a (N+l) trapped 
boson systems as [2(J Hl| 
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where Vt ra p(xi) is the external trapping potential and 
V{Si — Xj) is the two-body pair interaction. We 
use the standard Jacobi coordinates defined as Q = 

Xi+i - 7 J2]=i Xj , (i = 1) 2, ....N) and the cen- 
ter of mass through R — jyVy X)i=~i 1 x% ■ Then the relative 
motion of the atoms is described in terms of TV Jacobi 
vectors (Ci, • • •, Ov) as [H HI 



2i 
i+1 



pansion basis of the many-body wave function is the hy- 
perspherical harmonics (HH) . As HH basis keeps all pos- 
sible correlations, its direct application to trapped bosons 
in the condensate which contains at least few thousand 
bosons, is an impossible task. Very recently we have 
adopted a technique called potential harmonic expansion 
method (PHEM) which keeps all possible two-body cor- 
relations together with realistic interatomic interaction 
[2(|[22|. For the dilute Bose gas, the effect of two-body 
correlations are important and one can safely ignore the 
effect of all higher-body correlations. That is when two 
atoms interact the rest bosons are inert spectators and 
for zero temperature BEC this is a justified approxima- 
tion. Thus for the spinless bosons, we decompose "J in 
two-body Faddecv components 



N+l 



* = E ^ 



(3) 



ij>i 



Hence i/'ij is a function of two-body separation vector 
only, besides the global length (hyperradius, see below). 
ipij (symmetric under P^) satisfy the Schrodinger equa- 
tion 
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T being the total kinetic energy; operating JV on 
both sides of Eq. (4) we get back the original Schrodinger 

equation. The hyperradius is defined as r = \J^2f =1 Q ■ 
The hyperradius and (3TV — 1) hyperangles (denoted 
by fijv) together constitute 3TV hypcrspherical variables. 
The choice of Jacobi coordinates is not fixed as the la- 
beling of the particle indices is arbitrary. We choose a 
particular set for the (ij) interacting pair, called the (ij)- 
partition, by taking as £jv, and ($, tp) are two spherical 
polar angles of the separation vector fij . The angle <\> is 



defined through 



rcos< 



For the remaining (TV — 1) 



Jacobi coordinates we define the hyperradius for the par- 



tition (ij) as pij 
and 
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such that p: 



rsint 
coordinates are 



With this choice, the hyperspherical 



(r, £l N ) = (r, 4>, 1?, <p, f2jv-i ) 



(5) 



where f2jv— i involves (3TV — 4) variables: 2(TV — 1) polar 
angles associated with (TV— 1) Jacobi vectors £i, CiV-i 
and (TV — 2) angles defining the relative lengths of 
these Jacobi vectors [23]. Then the Laplacian in 3TV- 
dimensional space has the form 
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Hyperspherical harmonic expansion method (HHEM) is a 
convenient tool in many-body physics [23| , where the ex- 
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L 2 (f2jv) is the grand orbital operator in D = 3N dimen- 
sional space. Potential harmonics for the (ij)-partition 
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are defined as the eigenfunctions of L 2 (ilpf) correspond- 
ing to zero eigenvalue of L 2 (£l N _i) . The corresponding 
eigenvalue equation satisfied by L 2 (Qn) is [HI 

[L 2 (n N ) +C(C + D- 2)] T%g +l ((Uj) = 0, C = 2K+1- 

(7) 

This new basis is a subset of the HH basis and is called 
as potential harmonics (PH) basis. It docs not contain 
any function of the coordinate Q with i < N and is given 
by 



-yl,m 
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,(fi (tf) ) = Y lm {^ 3 ) {N) P l 2 Q K+l (</>)y (D - 3), (8) 

where Yi m (ujij) is the spherical harmonics and Wy = 

(&,<p)- The function W 'P^x+iifi) is expressed in terms 
of Jacobi polynomials and yo(3N — 3) is the HH of or- 
der zero in the (3iV — 3) dimensional space, spanned 

by {&}•'■ i0v-i} Jacobi vectors [23[. Thus the con- 
tribution to the grand orbital quantum number comes 
only from the interacting pair and the 3iV dimensional 
Schrodinger equation reduces effectively to a four dimen- 
sional equation. The relevant set of quantum numbers 
are three - orbital I, azimuthal m and grand orbital 2K+1 
for any N. The full set of quantum numbers are 

h =h = ■■■ = l N -i =0, l N = l (9) 
mi = m 2 = ■ ■ ■ = 97ijv-i = 0, mjy = m (10) 
n 2 = n 3 ~ ■ ■ ■ ~ ujv-i =0, n N ~ K. (11) 

We expand (ij)-Faddecv component, ipij, in the complete 
set of PH basis appropriate for the (ij) partition : 
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which includes only two-body correlations. Taking pro- 
jection of the Schrodinger equation on a particular PH, 
a set of coupled differential equations (CDE) is obtained 

MM 
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U K i(r) = fKiu l K (r), C = l+ 2 2 
I being the orbital angular momentum contributed by 
the interacting pair and K is the hyperangular momen- 
tum quantum number. Jki is a constant representing the 
overlap of the PH for interacting partition with the full 
set, which is given i n 12011 . The potential matrix element 
V K K' (r) is given by |20[ 



V KK >(r) 
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(Q^)n^)^«»- (14) 



So far we have disregarded the effect of the strong short 
range correlation in the PH basis. In the experimen- 
tally BEC, the average interparticle separation is much 



larger than the range of two-body interaction. This is 
indeed required to prevent atomic three-body collisions 
and formation of molecules. As the energy of the inter- 
acting pair is extremely small, the two-body interaction 
is generally represented by the s-wave scattering length 
(a s ). A positive value of a s gives a repulsive condensate 
and a negative value of a s gives an attractive conden- 
sate. However a realistic interaction, like the van der 
Waals potential, is always associated with an attractive 
— ?r tail at larger separation and a strong repulsion at 

ij 

short separation. Depending on the nature of these two 
parts, a s can be either positive or negative. For a given 
finite range two-body potential a s can be obtained by 
solving the zero-energy two-body Schrodinger equation 
for the wave function rj{xij) 
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+ V(x ij )r,(x ij ) = 0. (15) 



The correlation function quickly attains its asymptotic 
form C(l — ^-) for large Xij. The asymptotic normaliza- 
tion is chosen to make the wavefunction positive at large 
Xij (22)]. In the experimental BEC, the energy of the 
interacting pair is negligible compared with the depth 
of the interaction potential. Thus r)(xij) is a good ap- 
proximation of the short-range behavior of ipij . Then we 
introduce this correlation function in the expansion ba- 
sis and call it as correlated potential harmonics (CPH) 
basis. 
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r)(xij) correctly reproduces the short separation behavior 
of the interacting-pair Faddccv component. Convergence 
rate of the PH expansion is quite fast. The correlated 
potential matrix element Vkk> ( r ) is now given by 



V KK ,{r) = (hfhf,)- 1 * J + ^{pf(z)V^ 
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'\[^\ m(z)}dz. (17) 



Here h^f and wi(z) are respectively the norm and weight 
function of the Jacobi polynomial P%f(z) [2(| [HJ. K and 
K 1 are the grand orbital quantum numbers of the basis 
sets in which the potential matrix is calculated. 



III. CHOICE OF INTERACTION AND 
CALCULATION OF ENERGY LEVELS 

For the present study we consider few thousands 
(1000-10000) 87 Rb atoms in the JILA trap [H, [T3]. 

Throughout our calculation we choose a^ = J -|- as the 
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unit of length (oscillator unit) and energy is expressed in 
units of the oscillator energy (foui). The van der Waals po- 
tential has been chosen as the interatomic potential with 
a hard core of radius r c , viz, V(xij) = oo for Xij ^ r c 

and ~- for Xi 



> r c 



The strength C& is taken as 



6.4898 x 10" 11 o.u. for 87 Rb atoms in the JILA experi- 
ment [ItJ ■ We adjust the cut off radius r c in the two-body 
equation to correctly obtain a s = 2.09 x 10 -4 o.u. 

With these sets of parameters we solve the set of cou- 
pled differential equation [Eq. 13] by the hyperspherical 
adiabatic approximation (HAA) [24J. We assume that 
the hypcradial motion is slow in comparison with the hy- 
pcrangular motion. For the hyperangular motion for a 
fixed value of r, we diagonalize the potential matrix to- 
gether with the hypercentrifugal term. Thus for a fixed 
value of r, the equation for the hyperangular motion can 
be solved adiabatically. The eigen value of this equation 
is a parametric function of r and provides an effective 
potential for the hypcrradial motion. In the HAA, the 
lowest eigenpotential is used for the ground state of the 
system and the hyperangular motion appears through 
the coupling matrix Vkk' (t) . Thus the whole conden- 
sate collectively oscillates in the effective potential. The 
energy is thus obtained by solving the equation for the 
hyperradial motion as [24| 



m dr 2 



■w (r) 



E 

K=0 



dr 



E 



Co(0 = 0. 
(18) 

subject to the appropriate boundary conditions on 
Co(^*)' This is called uncoupled adiabatic approximation 
(UAA) . whereas disregarding the third term corresponds 
to extreme adiabatic approximation. HAA has already 
been successfully applied in different atomic and nuclear 
cases. 

Although the lower multipolarities have been 
successfully detected in the experiments [l8j|, however 
the collective excitations with higher multipolarity are 
also important specially to study the thermodynamic 
properties |19| . In our many-body picture, the collective 
motion of the condensate is characterized by the effective 
potential as described earlier. Thus the excited states in 
this potential are the states with the I th surface mode 
and nth radial excitation, which are denoted by E n i. 
Thus n=0 and 1=0 corresponds to the ground state and 
for I ^ 0, we get the surface modes. To calculate the 
higher levels with / ^Owc follow the next procdure. For 
I 0, a large inaccuracy is involved in the calculation of 
off-diagonal potential matrix and numerical computation 
becomes very slow. However, the main contribution 
to the potential matrix comes from the diagonal hy- 
percentrifugal term and we disregard the off-diagonal 
matrix element for I > 0. Thus we get the effective 
potential w;(r) in the hyperradial space for / ^ 0. The 
energy of the lowest modes are in close agreement with 
the other calculations fl8l . [l9l [25l . |26| and we observe 
that for energy much larger than the chemical potential, 



the excited states are separated at energy close to the 
harmonic energy tuo as in the noninteracting harmonic 
oscillator model. This transition from the low-energy 
collective to high-energy single particle spectrum leads 
us to study further the level fluctuation and other 
statistical behavior. 

Before discussing the statistical behavior of the 
energy spectrum we discuss how good is our approxima- 
tion method. There are many approximation methods 
to calculate the lowlying collective excitations and also 
the higher multipolarities. All these basically use the 
uncorrelated mean-field theory and the hydrodynamic 
method. Hydrodynamic method is good for large num- 
ber of bosons in the trap in Thomas Fermi limit. Where 
as our system is finite sized and have few thousands 
bosons. As the system is not exactly solvable like the 
ID system with contact S interaction, it is not possible 
to calculate the accuracy of our approximation method. 
However from our previous calculation of different 
measurable quantities like crtical instability of attractive 
BEC, the collective excitations, the thermodynamic 
properties we observed that the correlated PHEM is an 
improvement over the Gross Pitaevskii mean field treat- 
ment for several reasons. Although the GP mean-field 
equation is being widely used, the wave function does 
not include any correlation. It is pointed out by several 
authors that the replacement of the actual interaction 
by a contact potential is not appropriate for a general 
realistic potential which consists of a repulsive core 
and an attractive part [27M29I ] . The earlier studies also 
indicate the necessity of shape dependent potentials 
instead of the zero-range potential [10, HH . The choice 
of contact interaction specially for 3D attractive BEC 
is not satisfactory (30j as the ^-function interaction 
produces an essential singularity at r = 0. Thus to 
include correlations, one must go beyond the mean-field 
approximations and have to use the finite-range realistic 
potentials. Thus the correlated potential harmonics 
basis and the PHEM technique is a right step in this 
direction. This basis set retains all two-body correlations 
and assume that three and higher-body correlations 
are negligible. For dilute condensate it is perfectly 
justified. However including all two-body correlations 
we go beyond the mean-field theory. As the number 
of variables is reduced to only four for any number of 
bosons in the trap, we can treat quite a large number of 
atoms in the trap without much numerical complication. 
The use of van der Waals potential having a finite range 
takes care of the short-range repulsion and interatomic 
correlations. Clearly it is an improvement over the GP 
mean field theory. 

The correlated many-body aproach has been suc- 
cessfully applied in the calculation of static properties in 
different traps, collective excitations and thermodynamic 
properties of trapped bosons [HM H, H, H2IH. In 
Refs [2^, [26| , we have calculated the low- lying collective 
excitations by PHEM both for repulsive and attractive 
BEC. It has been shown that for repulsive and weak 
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interaction the many-body results arc close to the 
numerical solution of GP equation. However for large 
repulsive interaction significant difference is found. For 
the attractive BEC, the excitation frequencies for low 
lying modes are well comparable with the self-consistent 
Popov approximation. However higher multipolarities 
are lower than the GP result. This is attributed to the 
two-body correlations and finite-range interaction of the 
realistic interatomic interaction. In Ref [32|, the PHEM 



is extended to investigate thermodynamic quantities 
which involves the calculation of a large number of 
energy levels. The calculated critical temperature 
and the condensate fraction have been calculated and 
compared with the GP results. The effect of realistic 
interatomic interactions and two-body correlations 
on thermodynamic properties of trapped bosons are 
observed. Thus the calculated energy levels are accurate 
for further analysis. We check for the convergence such 
that the error is considerable smaller that the mean level 
spacing. 

Now to corroborate with the experiments we need 
the following discussions. Exciting the condensate by 
applying inhomogencous oscillatory force with tunable 
frequency, it is possible to observe several modes with 
different angular momentum and energy in the collective 
excitations [l8j . These experiments mainly concern 
low-lying collective excitations where the effect of 
interatomic interaction is prominant and the high-lying 
spectra should exhibit the single particle nature. The 
transition from collective to single particle excitations 
has also been studied in details theoretically The 
collective excitations have also dramatic dependence 
on the temperature which comes from the interaction 
between the condensate and the thermal cloud J3(| [I?} ■ 
But for the present study we consider the zero tem- 
perature BEC and the effect of thermal fluctuations 
do not arise. There is no damping in the condensate 
as we assume there are no thermally excited atoms. 
Apparently it may contradict the experimental situation. 
But in the presence of external trapping and at zero 
temperature the effect of damping is not critical. 



IV. STATISTICAL TOOLS AND NUMERICAL 
RESULTS 

After getting the many-body collective levels in- 
cluding higher order excitations with different Z, we trans- 
form the spectrum. Next to characterize the spectral 
fluctuation in the many-body energy spectrum and to 
compare the statistical properties of different parts of the 
spectrum we remove the smooth part in the level density. 
In general, the level density has two parts : one is the 
smooth part which defines the general trend of the en- 
ergy spectrum and a fluctuating part. The smooth part is 
removed by unfolding which maps the energy levels to an- 
other with the mean level density equal to 1 . Several un- 



folding procedures are in the Ref 0,111]. For the present 
calculation the many-body level density is approximated 
by a polynomial and unfolding is done by a 7th order 
polynomial. We unfold each spectrum separately for a 
specific value of I and then form an ensemble having 
the same symmetry. Next in order to study the spec- 
tral fluctuation of this unfolded spectrum we utilize the 
following statistical tools. Nearest neighbor spacing dis- 
tribution (NNSD) is the most applied tool in the study of 
the short range spectral correlations. From the unfolded 
spectrum, we calculate the nearest neighbor spacing as 
s = Ei + \ — Ei and calculate the probability distribution 
P(s). Uncorrelated spectrum obeys the Poisson statis- 
tics and P(s) = e~ s . Whereas for the system with time- 
reversal symmetry, level repulsion leads to Wigner-Dyson 

distribution P(s) = -^se^T - [39J. The A3 statistics is 
commonly used to investigate long-range correlations. It 
gives a statistical measure of the rigidity of a finite spec- 
tral level sequence. For a given interval L, it is often 
determined by the least square deviation of the staircase 
from the best straight line fits it. 

The P(s) distribution of the unfolded energy spec- 
trum of 1000 interacting bosons in the external trap is 
presented in Fig. 1 (a)-(d)[scc figure caption]. 
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FIG. 1: (Color online) The P(s) distribution is presented for 
different number of levels as indicated below the each panel, 
panel (a): lowest 100 levels, panel (b): 100 to 500 levels, 
panel (c): 500 to 1000 levels, panel (d): 1000 to 5000 levels. 
In each panel the histograms presents the P(s) distrubution 
for the Hamiltonian (1) with N = 1000 interacting bosons, 
the green dashed curve represents the Poisson distribution 
and the blue dotted curve represents the Wigner distribution. 
The magenta color solid line in the panel (d) corresponds 
to Brody distribution, corresponding Brody parameter being 
v = 0.04. 

For comparison Poisson statistics and GOE statistics 
are also given in the same figure. For lowest 100 levels 
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there is no level with very small spacing and no level 
beyond s = 2.0. Although the peak arises at s ~ 1.0, 
it strongly deviates from Wigner distribution. For such 
low- lying levels, the effect of interatomic interaction is 
dominating and levels are highly correlated. This is also 
intuitively right as we may write the many-body effective 
potential in the following way 



wo(r) =V frap + V(r). 

V(r) 



1 9 9 

>mcu r 



(19) 



Where Vt ra p is the external harmonic trap as described 
earlier and V(r) is obtained by the diagonalization of the 
potential matrix together with the hyper centrifugal re- 
pulsion. Now for small value of r (which corresponds to 
low- lying energy level) the effect of V(r) is dominating. 
Although it was expected that these levels would exhibit 
chaotic signature and follow Wigner distribution, but the 
level repulsion is masked due to the existence of external 
harmonic trap. Thus it exhibits a mixed statistics which 
could not be perfectly interpolated between Poisson and 
Wigner distribution by using Brody parameter H . Thus 
the evolution of P(s) distribution clearly shows the pres- 
ence of two energy scales. The situation becomes more 
interesting for intermediate levels. The P(s) distribu- 
tion for 100 to 500 levels exhibits two peaks as shown in 
Fig. 1(b). The first narrow peak appears at s = with a 
second broad peak near s = 1.5. For such intermediate 
levels, a part of the levels are correlated and shows nor- 
mal level repulsion, whereas the other part do not repel 
each other and try to maintain Poisson statistics which 
is reflected as a first peak at s = 0. The effect of level 
repulsion is manifested in the second peak. It is very 
similar to the classical mixed system, a part of phase 
space is completely regular with the other part chaotic. 
We have checked that by varying the number of levels 
in such an intermediate band, the width and peak val- 
ues change but qualitative features remain the same. For 
much higher levels, the effect of interatomic interaction 
gradually decreases and the effect of harmonic trap starts 
to dominate. Thus more and more states are coupled in 
regular uncorrelated distribution. It is well reflected in 
Fig. 1(c), where we see that a large part of levels try 
to exhibit Poisson statistics whereas a small fraction of 
levels is associated with a level repulsion, with a small 
peak near s = 2.0. For much higher levels the effect of 
interatomic interaction is almost negligible and the levels 
become regular and close to the integrable system. P(s) 
distribution is very close to Poisson, but the peak value 
at s = is less than 1. We fit the histogram with the 
Brody distribution [9], 



P(v, s) = (1 + v)as v exp{-as 1+u ) 



(20) 



where 



F 



( 2 + u 



1+v 



is the Brody parameter. 

The interesting feature of this distribution is that it inter- 
polates between Poisson distribution (y = 0) of regular 
systems and the Wigner distribution with v = 1. Thus 
the degree of chaos is determined by the value of v. For 



quantitative comparison, we fit the P(s) histograms to 
P(v,s) in Fig. 1(d) and the calculated repulsion param- 
eter is v = 0.04. 
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FIG. 2: (Color online) The P(s) distribution is presented for 
different number of levels as indicated below the each panel, 
panel (a): lowest 100 levels, panel (b): 100 to 500 levels, 
panel (c): 500 to 1000 levels, panel (d): 1000 to 5000 levels. 
In each panel the histograms presents the P(s) distrubution 
for the Hamiltonian (1) with iV = 5000 interacting bosons, 
the green dashed curve represents the Poisson distribution 
and the blue dotted curve represents the Wigner distribution. 
The magenta color solid line in the panel (d) corresponds 
to Brody distribution, corresponding Brody parameter being 
v = 0.05. 

The results for 5000 bosons is presented in Fig. 2(a)- 
(d)[sec figure caption]. As the condensate is repulsive, 
with increase in particle number, the condensate wave 
function spreads out as the net effective repulsion Na sc 
increases. With increase in interaction more and more 
many-body levels show level repulsion and we expect 
level spacing distribution close to Wigner which is 
very similar to the completely chaotic system. But in 
our system as the level repulsion is suppressed by the 
external trap, the P(s) distribution deviates from the 
Wigner distribution. This clearly shows the presence of 
two energy scales even for such intermediate levels. For 
500 to 1000 levels, we see quantum chaos sets in and 
P(s) is very close to Wigner distribution. To observe 
and determine the best fit window to the Wigner, the 
P(s) distribution for 501-600, 601-700 and 701-800 levels 
are plotted in Fig. 3. We observe that 601-700 is the 
best fit window and the corresponding Brody parameter 
is v = 0.9. However this energy window strongly 
depends on the number of atoms and the scattering 
length, whereas for much higher levels, we observe a 
crossover from Wigner like level repulsion to Poisson. 
In Fig. 2(d) we again compare the histogram with 
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FIG. 3: (Color online) The P(s) distribution is presented for 
different region of the spectrum, panel (a): 501-600 levels, 
panel (b): 601-700 levels, panel (c): 701-800 levels. In each 
panel the histograms presents the P(s) distrubution for the 
Hamiltonian (1) with N = 5000 interacting bosons, the green 
dashed curve represents the Poisson distribution, and the blue 
dotted curve represents the Wigner distribution. The black 
dot-dashed line corresponds to Brody distribution with the 



the Brody distribution and the corresponding Brody 
parameter is v = 0.05. 

We observe that P(s) distribution depends 
crucially on the number of levels and on the net effective 
interaction Na sc . To get more detail physical picture 
we calculate energy levels for 10000 bosons and plot 
the P(s) distribution in Fig. 4. Due to strong repulsive 
interaction (Na sc ) ~ 2.09, the low- lying levels are highly 
correlated and should show strong level repulsion. In 
Fig. 4, we present P(s) distribution for the lowest 50 
levels. Although the distribution has a sharp peak near s 
= 0.8, the distribution shows Wigner-like level repulsion. 
Comparing the same for 5000 bosons [Fig. 3], we observe 
the signature of level correlation for much lower levels. 
We also observe the uncorrelated Poisson distribution 
for higher levels as observed earlier. 



1 . 6 
1 . 4 
1.2 

1 

. 8 
. 6 
0.4 
0.2 





ft 



0.5 



1.5 

s 



2 2.5 3 3.5 



FIG. 4: (Color online) Plot of the P(s) distribution for lowest 
50 levels with 10,000 interacting bosons. Histograms repre- 
sent the many-body result obtained for the Hamiltonian (1). 
Green dashed curve represents the Wigner distribution. 

Now in this connection it worths to calculate 
the spectral distribution for the energy levels which are 
calculated by Gross-Pitaveskii mean-field equation. As 
the mean-field equation uses the zero-range contact in- 
teraction and it completely ignores the interatomic corre- 
lation, it is interesting to observe the effect of finite range 
interaction and intcrparticle correlation in the spectral 
statistics. For the calculation of energy levels we use the 
dispersion law of the discretized normal modes for spher- 
ical trap which is given by 



u>(n r , I) = uj(2n r 2 + 2n r l + 3n r + I) 2 . 



(21) 



n r is the radial quantum number and n r =0 corresponds 
to surface excitation, whereas the monopole oscillation 
corresponds to n r =l and 1=0. Note that Eq. (21) has 
the dependence only on the radial nodes and angular 
momentum, azimuthal degeneracy thus exists. Eq.(21) 
is also valid in the collisionless hydrodynamics where the 
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FIG. 5: (Color online) Fig. 5(a): Plot of the ground state 
energy (in o.u.) per atom (E/N), obtained from the many- 
body theory and GP equation, as a function of n r for I = 0. 
Fig. 5(b): The P(s) distribution obtained from GP calcula- 
tion is presented as Histogram. 



number of atoms in the trap is quite high. For our chosen 
set. N = 5000 and a sc = 2.09 x 10 -4 o.u., the parameter 
JVa " ~ 1.045, which is just greater than one and Eq.(21) 
be valid for lower n r and I. To see the accuracy of 
prediction [Eq.(21)], we plot in Fig. 5(a) the ground state 
energy per particle (E/N) as a function of n r for 1=0. 
The GP results start to be lower and lower for larger n r . 
The trend is maintained for other higher values of I. Our 
many-body results start to be higher near n r = 200 due 
to the presence of hyper centrifugal repulsion term in the 
many body equation [Eq. (13)]. As pointed in the Ref 
[40] that Eq. (21) is accurate for thj < p,. Thus for the 
calculation of P(s) distribution using GP results we take 
lowest 500 levels for which the above condition is valid. 
It guarantees that our choice of levels will be reliable for 
the calculation of spectral distribution. In Fig. 5(b) we 
plot the P(s) distribution for lowest 500 levels obtained 
from the dispersion law. It nicely shows the existence of 
large number of quasi-degenerate states as P(s) exhibits 
a sharp peak near s = 1. The existence of degeneracy 
is also seen in Eq.(21) where the discrete eigen mode 
frequency w(n r , I) of the spatial variation of density, ob- 
tained in the context of the hydrodynamic model of the 
condensate at low temperature, is a function of the ra- 
dial quantum no n r and orbital quantum no / only and 



hence is degenerate with respect to the azimuthal quan- 
tum no m. The contact <5-potential in the GP equation 
can not lift this degeneracy. There is a 5-type peak at 
about s = 1. 

Such 5-type peak in the P(s) distribution is called as 
Shnirelman peak [U |42[. In the year 1993, Shnirelman 
showed that systems with time reversal symmetry should 
exhibit the Shnirelman peak in the P(s) distribution. 
This behavior is again expected from Eq. (19). In the 
mean field results V(r) is calculated only taking the con- 
tact 5 interaction and ignores the interatomic correlation 
completely. Thus it can not lift the degeneracy of the ex- 
ternal harmonic trap completely. Whereas in our many- 
body calculation the short-range hard sphere below the 
cutoff radius and the -^f tail in the interatomic inter- 
action takes care of the effect of both short-range and 
long-range correlation and gives actual physical picture. 
But P(s) contains no information about spatial corre- 
lation. A simple measure of spacing correlation is the 
correlation coefficient C defind as 



c = E^- 1 )^+i- 1 )/E( s «- 1 ) 2 



(22) 



For uncorrelated spectra C = 0. The calculated value of 
C for Fig. 5(b) is 1.0 which again signifies the existence 
of bulk quasi-degenerate states. 

Here we observe that P(s) distribution strongly 
depends on the number of energy levels and also on 
the net effective interatomic interaction. Thus it is 
hard to say about correlation properties only from the 
study of P{s) distribution. It represents the study 
of the correlation properties in the large energy scale 
which will give new physical insight. The spectral 
rigidity is often used as a stronger tool than the level 
distribution in the analysis of complex systems as it can 
take into account of the long-range correlation between 
the levels while P(s) distribution takes into account 
only nearest neighbour correlations. So further studies 
are needed in this direction. We are mainly interested 
in A3 statistics of Dyson and Mehata [43| which gives 
a statistical measure of the rigidity of a finite spectral 
level sequence. For a level sequence with a constant 
average level spacing, the staircase function on the 
average follows a straight line. Thus A3- statistics gives 
a measure of the size of fluctuations of the staircase 
function about a best fit straight line. For Poisson 
spectrum, the levels are uncorrelated, the spectrum is 
rigid and (Aa(L)) = y|, whereas for GOE distribution, 
levels are strongly correlated and (A 3 (L)) oc log L. So 
to confirm our earlier findings in P(s) distribution we 
next study the correlation structure and A3 statistics 
for our system which are shown in Fig. 6 and Fig. 7. 
The results < A 3 (L) > are plotted against L for the 
same parameters and same number of levels as chosen 
in Fig. 1 and Fig. 2. Results for 1000 interacting bosons 
are presented in Fig. 6. For comparison we also plot 
< A 3 (L) > for Poisson and GOE. For low-energy levels 
we observe the same trend, bending towards the GOE 
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prediction. However the saturated value is well below 
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FIG. 6: (color onlne) Spectral average < As(L) > computed 
for the Hamiltonian (1) with 1000 interacting bosons in the 
external trap. Red color corresponds to lowest 100 levels. 
Green color corresponds to levels between 100 and 500. Blue 
color corresponds to levels between 500 and 1000. Violet color 
corresponds to levels between 1000 and 5000. The straight 
line corresponds to Poisson distribution and the black color 
corresponds to GOE result. 

the GOE prediction. It again reflects the fact that level 
repulsion due to interatomic interaction is screened due 
to the effect of external harmonic trap. For 1000-5000 
levels, < Aa(L) > follows the expected straight line 
behaviour upto L < 10, which is the result of integrable 
systems. But beyond L = 10, it still tends to saturation, 
but it is consistent with Berry's semiclassical arguments 
[idl |45| . The results for 5000 bosons are shown in the 
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FIG. 7: (Color online) Spectral average < Aa(L) > computed 
for the Hamiltonian (1) with 5000 interacting bosons in the 
external trap. Red color corresponds to lowest 100 levels. 
Green color corresponds to levels between 100 and 500. Blue 
color corresponds to levels between 500 and 1000. Violet color 
corresponds to levels between 1000 and 5000. The straight 
line corresponds to Poisson distribution and the black color 
corresponds to GOE result. 



Fig. 7, which again indicates that lower levels are highly 
correlated whereas for higher levels, we get the signature 
of more regular distribution. 

It is useful to mention that Guhr and Wei- 
denmuller [Irjj studied in the past the spectral properties 
of a regular Hamiltonain perturbed by a GOE. The 
results for 100 levels and 100 to 500 levels shown 
in the present work are quite similar to some of the 
results in Figs. 1,2,3 and 6 of [46| where a modified 
uniform spectrum was used as the regular Hamiltonian. 
Therefore, a quantitative description of the results in 
Figs. 1-5 in terms of a deformed GOE, which combines 
uniform, GOE and Poisson is possible but this is for a 
future investigation. This variety of behavior has also 
been observed early in the study of quantum mechanics 
of heavier clusters like Kr and Xe trimers. The energy 
spectrum of these clusters show a wide variety of 
behavior below and above the transition energy [47j . 
Very regular behavior of low-lying eigenstates changes 
to the combination of regular and irregular behavior at 
energy above the transition energy. 

Experiments on Bose-Einstcin condensation with 
7 Li atoms is another challenging research area where 
the s-wave scattering length (a s ) is negative, which 
indicates that atom-atom interaction is negative [48j |. 
A homogeneous condensate with a negative scattering 
length is impossible as the condensate approaches col- 
lapse. However the situation changes drastically in the 
presence of an external confinement. Spatially confined 
BEC is stable for a small, finite number of atoms (N cr ). 
For 7 Li, a s = -27.3 Bohr = -45.7xl0" 5 o.u. and for 
T = a metastable condensate exists when the number 
of atoms is less than the critical number ~ 1300 [49| . 
whereas theory predicts that BEC can occur in a trap 
with no more than about 1400 atoms [5(| ■ 

It is pointed out earlier in different connec- 
tion [5l|, [52| that the GP theory based on the pseudo 
potential form of the interatomic interaction is not 
suitable as an exact potential in 3D attractive system. 
Again as the attractive BEC becomes highly correlated 
near the critical points, the uncorrelated GP equation 
cannot take care of the effect of interatomic correlation. 
In our earlier calculation we have extensively applied our 
many-body method in the study of different properties 
and stability of the attractive condensate [53, [54|. In 
our present study we are interested in the spectral 
distribution of highly correlated BEC. The presence of 
hard sphere below some cutoff radius and the tail 
in the interatomic interaction properly takes care of the 
effect of both short-range and long-range correlations. 

For N < N cr , the condensate is metastable. In 
the many-body effective potential, the intermediate 
metastable region (MSR) is bounded by the high wall 
of the external trap on the right side and a very deep 
narrow attractive well appears on the left side of left 
intermediate barrier [54j . As the very high- lying levels 
will have large probability to tunnel through the barrier 
we are interested only for the low-lying levels for which 
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FIG. 8: (Color online) Plot of the P(s) distribution for dif- 
ferent number of levels in different region of spectrum with 
N — 1000 7 Li atoms. In each panel the histograms present 
the P(s) distribution obtained from many-body theory. The 
green dashed curve corresponds to the Wigner (GOE) dis- 
tribution in panel (a) and in panel (b) it corresponds to the 
Poisson distribution. 



the transition probability is almost zero. The results for 
N = 1000 and N = 1300 arc presented in Fig. 8 and in 
Fig. 9 respectively. Wc do not get any stable solution 
beyond TV = 1320. So for our present calculation the 
critical number is 7V cr =1320. For N = 1000, as the 
number of atoms is well below the critical number, we 
correctly obtain the lowest 100 levels which are well 
bound within the mctastable region. We observed high 
level correlation for lowest 100 levels and high-lying 
levels are uncorrelated. The effect of interatomic corre- 
lation for attractive BEC is very important for low-lying 
levels as the effect of V(r) dominates for smaller values 
of r. Although the level-correlation strongly dominates, 
however we do not observe any 6 function like peak. 
It signifies that for N = 1000, the exact degeneracy 
of harmonic oscillator is completely removed by the 
interatomic interaction V(r), whereas for larger r, as the 
term ^mui 2 r 2 dominates we get the uncorrelated Poisson 
distribution. The situation becomes more interesting 
for N = 1300 which is very close to the critical point 
and the condensate is highly correlated. It is reflected in 
Fig. 9(a)-(c) where we plot P(s) for different levels. As 
near the critical point, the metastable region becomes 
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FIG. 9: (Color online) The P(s) distribution for different 
number of levels of N — 1300 7 Li atoms in the trap is pre- 
sented as histograms. Panel (a) corresponds to lowest 50 lev- 
els, panel (b) corresponds to lowest 50-100 levels and panel 
(c) corresponds to lowest 100-200 levels as indicated in each 
panel. 



flatter, we calculate lowest 200 levels for which the 
tunnelling probability through the intermediate barrier 
is negligible. For lowest 50 levels [Fig. 9(a)] we observe 
a sharp peak in the first bin near s = 0. It signifies that 
many eigen-states overlap and it leads to the existence 
of large quasi-degenerate states. Although V(r) still 
dominates and the continuous symmetry of harmonic 
oscillator is removed, but discrete symmetry retains. In 
Fig. 9(b) and 9(c), although we get a sharp peak near 
s = 0, but it spreads to further bins. It signifies that 
the effect of quasi-degeneracy is gradually lifted in the 
higher levels. 
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V. CONCLUDING REMARKS 

Using the correlated many-body technique we com- 
pute the energy spectrum of weakly interacting trapped 
Bose gas. All throughout our calculation the Bose gas 
is at zero temperature and under the harmonic confine- 
ment with fixed trap size which corresponds to the JILA 
experiment. Although the statistical behaviour of com- 
pletely integrable and fully chaotic systems are under- 
stood, the intermediate region of integrability and chaos 
is more interesting. Interacting trapped bosons is such a 
system which is spatially inhomogcneous. The existence 
of external harmonic trap together with interatomic in- 
teraction makes the system more interesting. We study 
the spectral fluctuation and level correlation in the en- 
ergy spectrum. Although there is no rigorous derivation, 
but the numerical results show a mixed statistics, which 
is very complexly dependent on the number of energy 
levels and number of the bosons in the trap. However for 
higher energy levels where the external trap is dominat- 
ing we get back the Poisson type fluctuation, whereas the 
low-lying collective excitations are strongly influenced by 
interatomic interaction and shows level repulsion. Thus 
our findings do not strictly obey the earlier conjecture of 
Bohigas for atomic nucleus and atoms. The results for 



attractive Bose gas near the critical point is highly in- 
teresting which nicely shows how the degeneracy of the 
harmonic trap is gradually removed for higher levels. Al- 
though there is no experimental data for such high-lying 
states, but for dilute interacting Bose gas it is possible 
to measure experimentally with present-day set up. Our 
present study opens many questions for further study. In 
the present study the interacting Bos gas is in a fixed trap 
size. However the use of time-dependent potential will al- 
low to study the dynamical behavior of energy spectrum 
and the time evolution of the spectral statistics and corre- 
lation properties. Throughout our calculation we use the 
zero-temperature Bose gas and delibarately avoids any 
thermal fluctuations. So it is also interesting to study 
the spectral distribution for non zero temperature BEC. 
Our present methodology is not valid as it avoids the 
thermal fluctuation. 
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